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Abstract. The main aim of this paper is to document the performance of p-refinement with 
respect to maximum principles and the non-negative constraint. The model problem is (steady- 
state) anisotropic diffusion with decay (which is a second-order elliptic partial differential equation). 
We considered the standard single-field formulation (which is based on the Galerkin formalism) and 
two least-squares-based mixed formulations. We have employed non-uniform Lagrange polynomials 
for altering the polynomial order in each element, and we have used p = 1, • • • ,10. It will be shown 
that the violation of the non-negative constraint will not vanish with p-refinement for anisotropic 
diffusion. We shall illustrate the performance of p-refinement using several representative problems. 
The intended outcome of the paper is twofold. Firstly, this study will caution the users of high-order 
approximations about its performance with respect to maximum principles and the non-negative 
constraint. Secondly, this study will help researchers to develop new methodologies for enforcing 
maximum principles and the non-negative constraint under high-order approximations. 



1. INTRODUCTION 

Contaminant transport, chemical and biological remediation, and carbon-dioxide sequestration 
are some of the main engineering challenges of the 21st century. A large-scale implementation 
of any of these problems will have irreversible consequences on the environment, and can affect 
large geographical region and for long periods of time. In addressing these pressing issues, robust 
predictive numerical simulations have played and will continue to play an important role. 

In these kind of problems, predicting the fate of chemical species is an important component, 
and diffusion is a dominant phenomenon. It should be noted that concentration (and certain 
other quantities such as absolute temperature, density, absolute pressure, saturation in multiphase 
systems) attains only non-negative values. A negative value for the concentration is unphysical. In 
a coupled reactive-transport numerical simulator, a negative value for the concentration of a species 
will result in an algorithmic failure. A robust numerical simulation should therefore preserve the 
physical and mathematical requirement of non-negativeness of species concentration. 
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1.1. Maximum principles and high-order approximations. Mathematically speaking, steady- 
state diffusion- type equations are elliptic partial differential equations, and are known to satisfy 
the so-called maximum principles [TO]. The non- negative constraint can be obtained from maxi- 
mum principles under certain assumptions on the input data. The discrete version of maximum 
principles are commonly referred to as discrete maximum principles (DMP). A study on discrete 
maximum principles investigates whether a given numerical formulation (that is, in a discrete set- 
ting) inherits the underlying maximum principles (which are satisfied in the continuous setting). 
A robust predictive numerical simulation should preserve various fundamental properties like the 
non-negativeness and maximum principles. Hence, one is particularly interested in the necessary 
and/or sufficient conditions under which a given numerical formulation satisfies discrete maximum 
principles and the non- negative constraint. 

Many numerical formulations (under finite element, finite volume, and finite difference method- 
ologies) have been developed. However, it should not be expected that any of these formulations will 
satisfy maximum principles and the non-negative constraint as there are no in-built mechanisms 
in these formulations to meet such constraints. Recently, researchers have proposed numerical 
methodologies for enforcing maximum principles and the non-negative constraint. For example, 
References [19} [20] have addressed the non-negative constraint under the finite volume method. 
Liska and Shashkov |21] have proposed a methodology to meet these constraints using the conser- 
vative finite difference technique. In References [231 122] . optimization-based techniques have been 
employed to meet these constraints under both mixed and single-field finite element formulations 
but have restricted their studies to low-order finite elements. 

One of the early studies on discrete maximum principles is due to Varga [34J , which was in the 
context of finite difference. An important work on discrete maximum principles with respect to 
the finite element method is by Ciarlet and Raviart [9]. In this paper, the authors have shown 
that acute-angle triangulation is a sufficient condition for a low-order approximation to satisfy 
maximum principles and the non-negative constraint. However, this condition is not sufficient under 
a high-order approximation. Other notable works on discrete maximum principles for high-order 
approximation are [14. 361 l3T] I30 [ 135] . All these works considered one-dimensional problems. Herein 
we systematically study the performance of high-order approximations on several two-dimensional 
problems. 

Also, earlier numerical works on discrete maximum principles concentrated on single-field formu- 
lations [22], and mixed formulations based on the variational multiscale formalism or lowest-order 
Raviart-Thomas spaces [23]. Herein, we shall also investigate the performance of mixed formula- 
tions based on the least-squares formalism. 
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1.2. Least-squares formulations. Least-squares variational methods often constitute an appeal- 
ing alternative to the more popular weak formulation based on the Galerkin formalism in develop- 
ing efficient and robust finite element models. In a least-squares-based formulation a non-physical 
least-squares functional is defined in terms of the sum of the squares of appropriate norms of the 
governing partial differential equation residuals. Since the least-squares functional is by definition 
positive and convex and it naturally follows that the minimizer of the least-squares functional 
coincides with the exact solution of the original set of partial differential equations. The finite 
element model associated with the least-squares formulation is constructed via direct minimization 
of the least-squares functional with respect to the trial space associated with the finite element 
discretization. Some representative references on least-squares formulations are [H [16] . 

1.3. Main contributions of this paper. The main aim of this paper is to document the per- 
formance of p-refinement for solving diffusion-type equations with respect to maximum principles 
and the non-negative constraint. We illustrate the performance on various computational grids 
and using different canonical problems. We consider three weak formulations: the standard single- 
field formulation and two least-squares-based formulations. We illustrate the extent to which these 
formulations violate maximum principles and the non-negative constraint under high-order approx- 
imations. We do not, however, provide a methodology for enforcing maximum principles and the 
non- negative constraint under high-order approximations. We also show that the performance of 
least-squares formulations with respect to the non-negative constraint depends on the choice of 
the weight in the inner product. The present work has two main purposes. First, users of high- 
approximations will be aware of their performance with respect to maximum principles and the 
non-negative constraint. Second, it will help researchers to develop methodologies for enforcing 
maximum principles and the non- negative under high-order approximations. 

1.4. An outline of the paper. The remainder of this paper is organized as follows. In Section 
[2j we present the governing equations of anisotropic diffusion with decay. We shall also discuss 
the classical maximum principle, and its consequences (in particular, the non-negative constraint). 
In Section [3j we present the weak formulations, and in Section [4] we discuss high-order spectral 
approximations. In Section [5j we shall show the performance of these numerical formulations using 
several canonical problems, and conclusions are drawn in Section [6} 

We shall employ the following symbolic notation in this paper. To distinguish vectors in the 
continuum setting from vectors in the finite element context, we shall employ lower case boldface 
normal letters for the former, and lower case boldface italic letters for the later. For example, x is 
used to denote a spatial position vector, and c is used to denote a finite element vector containing 
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nodal concentrations. To distinguish a second-order tensor from a matrix, we shall denote second- 
order continuum tensors using upper case boldface normal letters, and shall denote matrices using 
upper case boldface italic letters. For example, D is used to denote the diffusivity tensor, and K is 
used to denote stiffness matrix. Throughout this paper, repeated indices do not imply summation. 
(That is, Einstein's summation convention is not employed.) We shall denote the set of natural 
numbers as N, and the set of real numbers as R. Other notational conventions are introduced as 
needed. 

2. GOVERNING EQUATIONS: DIFFUSION WITH DECAY 



In this paper, we assume that O is an open bounded subset of K , where "no?" denotes the 
number of spatial dimensions. The boundary is denoted by dVt := Cl — Q, where a superposed bar 
denotes the set closure. The boundary is divided into two parts: T D and T N , where T D is that 
part of the boundary on which Dirichlet boundary conditions are prescribed, and T N is the part 
of the boundary on which Neumann boundary conditions are prescribed. For well-posedness, we 
have r D U T N = and T D n T N = 0. A spatial point is denoted by x £ H. The gradient and 
divergence operators with respect to x are, respectively, denoted by grad[-] and div[-]. Herein, we 
shall consider the anisotropic diffusion of a chemical species, and allow the decay of the chemical 
species. Let us denote the concentration of the chemical species by c(x). We present the governing 
equations in the standard divergence form, which take the following form: 

a(x)c(x) — div[D(x)grad[c]] = /(x) in £2 (la) 
c(x)=c p (x) onr D (lb) 

n(x) • D(x)grad[c] = t p (x) on T N (lc) 

where a(x) > is the decay coefficient, D(x) denotes the diffusivity tensor, /(x) is the volumetric 
source/sink, c p (x) is the prescribed concentration, i p (x) is the prescribed flux, and n(x) is the unit 



outward normal vector on the boundary. The mathematical model given by equations (la)-(lc 
frequently arises in Mathematical Physics (see the discussion in |22l Introduction]). 

We shall assume that the decay coefficient is bounded above, which means that there exists a 
real constant ceo < +°o such that we have 

a(x) <a Vx 6 ft (2) 

(Note that we have already assumed that the decay coefficient is non- negative.) The diffusivity 
tensor is assumed to be symmetric, bounded above and uniformly elliptic. That is, there exists two 
constants < Ai < A2 < +00 such that 

Aiy T y < y T D(x)y < A 2 y T y Vx £ U, Vy G R nd (3) 
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It should be noted that uniform ellipticity condition is a stronger requirement than demanding 
that the matrix corresponding to tensor D(x) being positive definite at every spatial point in the 



domain. The system of equations (la)-(lc) is a second-order elliptic partial differential equation 
with Dirichlet and Neumann boundary conditions. It is well-known that the above partial differ- 
ential equation satisfies the so-called maximum principles under certain regularity assumptions on 
the domain and on the input data (D(x), a(x), /(x), c p (x), and t p (x)). 

2.1. Maximum principles and the non-negative constraint. There are different kinds of 
maximum principles for second-order elliptic partial differential equations available in the literature 
(for example, see the discussion in Han and Lin [12]). Herein, we consider two maximum principles 
by E. Hopf |15j . (Also see the commentary on Hopf's paper by J. Serrin [29j.) The first maximum 
principle is for pure diffusion without decay (that is, a(x) = 0), and the second maximum principle 
allows the possibility of decay (that is, a(x) > 0). For k £ N, we shall use C fc (f2) to denote the 
set of functions having all derivatives up to the order k continuous on Q, and C (O) to denote the 
space of all continuous functions on O that can be continuously extended to the boundary. 

Theorem 1 (A maximum principle for diffusion equation). Let c(x) £ C 2 (il) n C°(0) satisfy the 
following differential inequality 

-div[D(x)grad[c]] = /(x) < 
with D(x) is uniformly elliptic, bounded above and continuously differentiate. Then, we have 

max c(x) < max c(x) 

Theorem 2 (A maximum principle for diffusion with decay). Let c(x) £ C 2 {VL) n C°(Cl) satisfy 
the following differential inequality 

a(x)c(x) - div[D(x)grad[c]] = /(x) < 

with a(x) > 0; D(x) is uniformly elliptic, bounded above and continuously differentiable; and a(x) 
is bounded. Then, we have 

max c(x) < max {c(x),0} 

Mathematical proofs to the above two theorems can be found in Gilbarg and Trudinger [10]. In 
both theorems, we have assumed volumetric sink (i.e., /(x) < 0). If we have volumetric source 
(i.e., /(x) > 0) then the "max" has to be replaced by "min." Put differently, for volumetric source, 
the minimum occurs on the boundary in the case of pure diffusion, and the non-negative occurs on 
the boundary in the case of diffusion with decay. Similarly, if /(x) = 0, then both the maximum 

and minimum occur on the boundary in the case of pure diffusion, and the non-negative maximum 
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and non-negative minimum occur on the boundary in the case of diffusion with decay. If there are 
regions of both source and sink in the domain then the above theorems do not directly apply. The 
discrete version of maximum principles is commonly referred to as discrete maximum principles. 

There are several studies on discrete maximum principles, and some representative ones are 
[SI Q21 Q31 [13 HU E3[ E2] • Some of these studies (particularly the ones that derive necessary and 
sufficient conditions to meet the non-negative constraint and maximum principles, for example, 
references [18]) hinge on the following mathematical arguments. After spatial discretization using 
the finite element method (or, for that matter, using the finite volume method or finite difference 
method) one obtains a system of linear equations of the following form: 

Ax = b (4) 

where A is called the coefficient matrix, and b is called the force matrix. Under the low-order 
approximation (i.e., p = 1), it can be shown that b y if the forcing function /(x) > 0. (Herein, 
we have used the symbol y to denote component- wise inequality. That is, b y means that bi > 0.) 
If b y 0, a sufficient condition for x y is requiring the matrix A to be monotone (which means 
that the inverse of A has all non- negative entries). Usually, a stronger condition requiring that 
the matrix A is an M- matrix is commonly employed [18J. (There are different ways of defining an 
M- matrix, and for further details see references [21 EH]-) Restrictions are then placed (typically on 
the mesh) to make A to be an M- matrix. However, as discussed in references [23] and [22]) such 
a procedure will work only for isotropic diffusion, and is not sufficient for anisotropic diffusion. 
Several examples are presented in these two references to illustrate this conclusion. 

For high-order approximations, another complexity arises due to the fact that even if /(x) > 
Vx G its L2 projection onto the polynomial space need not be non-negative. That is, even if 
/(x) > 0, the condition 6^0 need not hold under high-order approximations. Hence, the whole 
analysis requiring that A to be a monotone or an M-matrix will not guarantee that the non- negative 
constraint and maximum principles under high-order approximations will be met. In this paper, we 
carefully assess the performance of high-order approximations with respect to maximum principles 
and the non- negative constraint. 

Remark 3. Some test problems in Section [5| do not have classical solutions in the sense that 
c(x) E C 2 (f2)nC°(il), and hence, theorems^ and^ do not directly apply. In the literature, however, 
one can find maximum principles for weak solutions posed under weaker regularity assumptions. For 
example, see references [251 [TJ [5] . Such a detailed discussion is beyond the scope of this paper, and 
is not central to our presentation. 

Remark 4. Some prior numerical works have addressed the non-negative constraint and maxi- 
mum principles for models similar to equation ([!]) on general computational grids. Nakshatrala and 
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Valocchi [23J have considered discrete maximum principles for pure diffusion and two mixed for- 
mulations. Specifically, they considered the variational multiscale formulation and the lowest-order 
Raviart- Thomas element. They provided a methodology for enforcing maximum principles and the 
non-negative constraint that works only for low-order finite elements. Nagarajan and Nakshatrala 
|22j have considered discrete maximum principles for the diffusion equation with decay under the 
standard single-field formulation. They also restricted their studies to low-order finite element 
approximations. Liska and Shashkov [2Tj have addressed discrete maximum principles for pure dif- 
fusion using conservative finite difference techniques. The present study differs from previous works 
in two ways: use and performance study of (I) high-order approximations, and (2) least-squares 
finite element models. 

3. WEAK FORMULATIONS 

In this paper we shall consider three weak formulations. The first is the classical single-field 
formulation, which is based on the Galerkin formalism. The second and third are mixed least- 
squares-based formulations. For completeness, we briefly outline these formulations. To this end, 
let us introduce the following function spaces: 

U := (c(x) € H X (VL) I c(x) = c p (x) on T D } (5a) 
W:= {w(x)eff 1 (n)|w(x)=Oonr D } (5b) 

Q:={q(x)|q(x)G(if 1 (^)r d } (5c) 

where H 1 ^) is a standard Sobolev space [32]. (Recall that "nci" denotes the number of spatial 
dimensions.) 



3.1. The classical single-field formulation. The weak formulation based on the Galerkin for- 



malism for the governing equations (la)-(lc) can be stated as follows: Find c(x) 6 IA such that we 
have 



B(w(x); c(x)) = T(w(x)) Vuj(x) € W 



(6) 



where the bilinear form B(w(x); c(x)) and the linear functional J-(w(x)) are, respectively, defined 
as follows: 

B(w(x); c(x)) := / a(x)u>(x)c(x) dfi + / grad[w] • D(x)grad[c] dSl (7a) 
Jn Jn 

.F(w(x)) := / w(x)/(x) dfl+ [ w{x)t p {x) dT (7b) 



One can easily verify that the above weak formulation ( 23 ) is equivalent to minimizing the following 
functional: 

J(c(x)) := i J (a(x)c 2 (x) + grad[c] • D(x)grad[c])dft - f c(x)/(x) dft 

-[ c(x)tP(x)dr=^( C (x); C (x))-^(c(x)) (8) 

Hence, the finite element approximation inherits the best approximation property with respect to 
the norm ||c(x)||£ = J |£>(c(x); c(x)). 

3.2. Least-squares-based formulations. Least-squares variational procedures possess many at- 
tractive mathematical as well as practical computational attributes. In particular, the least-squares 
method always invokes a minimization principle whose minimizer coincides with the solution of the 
governing partial differential equations. As a result, the discrete numerical solution always pos- 
sesses the best approximation property with respect to a well-defined norm (i.e., the energy norm 
of the functional). When this norm is equivalent to a standard norm (such as a norm from an 
appropriate Sobolev space) ideal convergence rates of the finite element solution may be estab- 
lished. That the least-squares formulation is always based on a minimization principle insures a 
robust setting that is often lacking in the associated weak form Galerkin model. In addition, the 
bilinear form resulting from invocation of the minimization principle in least-squares formulations 
will always be symmetric and positive definite, and restrictive compatibility conditions such as the 
inf-sup condition will never arise. As a result, providing a conforming discretization and well-posed 
boundary conditions, the discrete variational formulation will always possesses a unique solution. 
Furthermore, the symmetry and positive-definiteness of the resulting coefficient matrix may be ex- 
ploited directly in the solution process of the linear system of equations. In particular, direct solvers 
may employ a sparse Cholesky decomposition of the coefficient matrix, while iterative solvers may 
utilize the robust preconditioned conjugate gradient algorithm. Due to symmetry, only half of the 
coefficient matrix need to be actually calculated and stored in memory. 

Least-squares formulations are not without their own deficiencies and in many cases there are 
fixes. For example, unlike weak form Galerkin formulations where regularity requirements of the 
finite element spaces are weakened (by invoking Green's identities), least-squares formulations re- 
quire higher regularity of the finite element spaces (dictated by the order of the governing partial 
differential equations). Higher regularity requirements negatively affect the condition number of 
the coefficient matrix and also the continuity requirement of the solution across element bound- 
aries. High regularity requirements may be avoided by constructing the least-squares finite element 
model in terms of an equivalent lower-order system by the introduction of additional independent 
variables. Such a mixed formulation allows for the use of standard Lagrange interpolation functions 



but also produces a large set of global equations to be solved. Such a formulation is often quite 
valuable, however, since the auxiliary variables are typically physical quantities of interest such as 
the flux. 

A least-squares-based formulation will be based on the minimization of a least-squares functional, 
which is constructed from the sum of the squares of appropriate norms of the residuals of the 
governing partial differential equations and the Neumann boundary condition. Although it is 



certainly possible to construct a least-squares finite element model based on equations (la)-(lc 
directly, such an approach requires a higher degree of regularity in the finite element solution such 
as c(x) £ H 2 (fi), where H 2 (Q) is another standard Sobolev space on |32j . Herein, we recast the 
governing equations into the following equivalent first-order system: 

a(x)c(x) + div[q] = /(x) in O (9a) 

q(x) = -D(x)grad[c] in O (9b) 

c(x) = c p (x) onr D (9c) 

-n(x) • q(x) = i p (x) onr N (9d) 



where q(x) is the diffusive flux vector field. We shall treat both c(x) and q(x) as independent 
variables, and employ the standard Lagrange finite element interpolation functions for both of 
these field variables. The price incurred will be an increase in the overall size of the system of 
equations that one needs to solve to obtain a numerical solution. The additional expense is not 
unwarranted as the flux is an important quantity of interest in many engineering applications. As a 
result, one obtains the flux directly in the solution process as opposed to the usual way of obtaining 
it at the post-processing stage of a finite element analysis. 

Remark 5. For stability reasons, some least-squares formulations for diffusion-type equations aug- 



ment the first order form of the governing partial differential equations (i.e., equations (9a)-(9d)j 
with the following seemingly redundant expression: 

curlfq] = (10) 

For example, see Reference [4J. It is important to note that q satisfies the above expression, in 
general, only if the medium is homogeneous ( that is, the diffusivity tensor is independent ofx) and 
isotropic. In this work, we place no such restrictions on D(x), and hence do not utilize the above 
expression in construction of the least-squares functional. 
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We shall define the norm II • \\k for scalar and vector functions defined over K as follows: 



x)[&:= / u 2 (x)dK (11a) 

JK 



IN 

\\u(x)f K := / u(x)-u(x)dK (lib) 
Jit 

In this paper we investigate the performance of two mixed least-squares-based formulations, and 
denote them as "LSI formulation" and "LS2 formulation". For both these formulations, the least- 
squares functional can be constructed as follows: 

JLS (c(x), q(x)) := \ ||/3(x) (a(x)c(x) + div [q] - /(x))|| 2 



1 



+ - || A(x) (q(x) + D(x)grad[c] 



In 



+ ^||q(x)-n(x) + iP(x)||2 N (12) 
where the second-order tensor A(x) and the scalar function /3(x) are, respectively, defined as follows: 



[ I LSI formulation 

A(x) = { , (13a) 
D _1 / 2 f x ) LS2 formulation 



/3(x) 



1 LSI formulation 

a _1 / 2 (x) ifa(x)^0 



1 if a(x) = I (13b) 
> LS2 formulation 



where I is the second-order identity tensor. Recall that a(x) > 0, and hence /3(x) is well-defined. 
Also, note that D(x) is a positive definite tensor, and by square root theorem D -1//2 (x) is also 
well-defined |llj . The variational principle associated with both the least-squares formulations can 
be stated as follows: Find c(x) £ U and q(x) G Q such that we have 

J hS (c(x), q(x)) < Jls (c(x), q(x)) Vc(x) G ZY, q(x) G Q (14) 

The first-order optimality condition demands that the first variation of ^ls ( c ( x )) f l( x )) be- iden- 
tically zero. The corresponding weak statement can be written as follows: Find c(x) G tl and 
q(x) G Q such that we have 



B LS (ifl(x),w(x);c(x),q(x)) = T ts (to(x), w(x)) Vw(x) G W, w(x) G Q 
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(15) 



where the bilinear form £>ls (^( x )> w ( x ); c(x), q(x)) and the linear functional Jxs {w(x), w(x)) are, 
respectively, defined as follows: 

fi L sMx),w(x);c(x),q(x)) := / /? 2 (x) (a(x)w(x) + div[w]) (a(x)c(x) + div[q]) dft 

+ / (w(x) +D(x)gradH) ■ A 2 (x) (q(x) + D(x)grad[c]) dft 

+ / w(x) • n(x) q(x) • n(x) dr (16a) 

J"ls(w(x),w(x)) := / /3 2 (x) (a(x)w(x) +div[w])/(x) dfi - / w(x) • n(x) t p (x) dr 
Jn ' Jr N ' 

(16b) 

4. HIGH-ORDER APPROXIMATIONS 

As shown in the previous section, a variational formulation (either Galerkin or least-squares 
formulations) of a general boundary value problem may be stated as follows: Pind u(x) £ V such 
that we have 

B(w(x);u(x)) = J r (w(x)) Vw(x)GV (17) 

where i3 (w(x); u(x)) is a bilinear form, ,F(w(x)) is a linear form, and V and V are appropriate 
function spaces. The quantity u(x) represents the set of independent variables (associated with the 
variational boundary value problem), and w(x) represents the corresponding weighting function. 

In the finite element method we restrict the solution space to a finite dimensional sub-space V hp 
of the infinite dimensional space V, and the weighting functions to a finite dimensional sub-space 
yh%, c y As a result, in the discrete case we seek to find u^ p (x) £ V hp such that we have 

B(w ftp (x);u ftp (x))=^(w hp (x)) Vw hp (x) € V fcp (18) 

We assume that the domain Q C M. nd is discretized into a set of NE non-overlapping sub-domains 

Cl e , called finite elements, such that Q, ~ Q hp = U^=i The geometry of each VL e is characterized 

using the standard isoparametric bijective mapping from Cl e to the master element £l e . In the 

present study we restrict the classes of elements considered to lines in R 1 , four side quadrilateral 

elements in R 2 and six face brick elements in R 3 (although numerical results are presented for 

nd = 1 and 2 only). As a result we can simply define the master element as Cl e = [— l,+l] nd . 

The natural coordinates associated with f2 e (when nd = 3) are defined as £ = (£, n, £). In this 

work we utilize a family of finite elements constructed using high polynomial order interpolation 

functions. The quantities h and p appearing in the definition of the sub-space V hp imply that the 

discrete solution may be refined by either increasing the number of elements in Q hp (/i-refinement), 

increasing the polynomial order of the approximate solution within each element Cl e (p-refinement) 

or through an appropriate and systematic combination of both h- refinement and p- refinement. 

li 



Within a typical finite element a given variable, such as the species concentration c(x), may be 
approximated by the following formula 

n 

c (x) « c hp (x) = (0 in ^ (19) 

i=l 

where ipi (£) are the nd-dimensional Lagrange interpolation functions, cf are the values of Ch p (x.) at 
the element nodes and n = (p + l) nd is the number of nodes in Cl e . In addition, the quantity p is the 
polynomial order of each interpolation function. There are a variety of ways in which high-order 
nd- dimensional interpolation functions may be formulated. For our analysis we construct these 
polynomial functions from tensor products of the one-dimensional C° spectral nodal interpolation 
functions 

« ( « = P ( P + i)Mf 3 )tt-« m[ " 1 ' 11 (20) 

where L p (£) is the Legendre polynomial of order p. The quantities £j represent the locations of the 
nodes associated with the one-dimensional interpolants (with respect to the natural coordinate £). 
The one-dimensional nodal points are defined as the roots of the following expression 



(£-l)(£+ 1)^(6 = in [-1,1] (21) 

Figure [T] compares high-order interpolation functions for p = 6 using uniform nodes with corre- 
sponding interpolation functions using non-uniform nodes based on Gauss-Lobatto-Legendre points. 

Multi-dimensional interpolation functions may be constructed from simple tensor products of the 
one-dimensional spectral interpolants. For example, in two-dimensions, the high-order interpolation 
functions may be defined as 

M^v) = 'Pj(.0<Pk(v) in^ e = [-l,l] 2 (22) 

where i = j + (k — 1) (p + 1) and j,k = 1, • • ■ ,p + l. Likewise, in three-dimensions, the interpolants 
are expressed as 



& (£. V, = <Pj (0 <Pk (V) VI (0 m n e = [-1, l] 3 (23) 

where i = j + [k — 1 + (I — 1) (p + 1)] (p + 1) and j, k, I = 1, • • • ,p + 1. Finite elements constructed 
from tensor products of </?j (£) are commonly referred to as spectral elements in the literature [17] ■ 
Such elements are merely standard high order Lagrange type finite elements, where the locations 



of the unequally spaced nodes in Q e are taken as tensor products of the roots of equation (21). 

Figure [2] illustrates various examples of high-order two-dimensional elements. 
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The finite element formulation naturally leads to a set of linear algebraic equations. The resulting 
coefficient matrix and force vector associated with the equations are constructed from the bilinear 



and linear forms appearing in equation (18) respectively. In this work we utilize the standard 
Gauss-Legendre quadrature rules for numerical integration of these quantities. We utilize full 
integration of all integrals and do not resort to selective under-integration of any terms in the 
coefficient matrix or force vector. All numerical results have been obtained using a quadrature rule 
of at least NGP = p + 1, where NGP represents the number of quadrature points in the direction 
of a given natural coordinate associated with f2 e . For details on the computer implementation of 
the finite element method, including descriptions of the bijective isoparametric mapping Cl e «^ Cl e 
and the global assembly operator, we refer to the texts of Reddy [26] and Bathe [3]. For details 
on construction of the spectral interpolation functions, we refer to the book by Karniadakis and 
Sherwin (TT] . 

5. REPRESENTATIVE NUMERICAL RESULTS 

In this section, we illustrate the performance of high-order approximations under the classical 
single-field formulation and two least-squares-based formulations with respect to maximum princi- 
ples and the non-negative constraint. Before we present numerical results, we briefly discuss how 
the numerical results under high-order approximations are visualized. 

5.1. Visualization of results using high-order approximations. Once a numerical solution 
has been obtained using a particular finite element discretization, a given field variable may be 
evaluated at any point within a typical finite element using the standard interpolation formula 



given in equation (19). Use of this formula is crucial in evaluating the performance of high-order 
finite element models with respect to maximum principles and the non- negative constraint. Un- 
fortunately, it is typically impossible to employ this scheme explicitly in the actual visualization 
of multi-dimensional numerical results obtained using high-order spectral//ip finite element formu- 
lations. Visualization software (such as Tecplot PQ) allow for the post-processing of structured 
data associated with a given finite element mesh; however, such programs typically require data 
structures containing the element connectivity array of low-order elements only. To utilize standard 
visualization software, it therefore becomes necessary to convert data associated with a high-order 
spectral/fop finite element mesh into data associated with a low-order finite element mesh. 

Solution data associated with a high-order finite element mesh may be readily converted into 
solution data associated with a low-order finite element discretization through the creation of a 
set of fictitious low-order visualization elements, utilized for plotting purposes only. There are a 
variety of ways in which a low-order visualization mesh may be created. Perhaps, the simplest 
choice is to create a low-order mesh using only the actual nodes of the high-order finite element 
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discretization. Unfortunately, the visualized numerical results will inevitably deviate from the 
actual finite element solution within a given element when such an approach is taken. In an effort 
to minimize visualization errors in the presentation of our numerical results, we utilize equation 



( 19 ) to evaluate the numerical solution at a discrete number of points within a given element that 
is greater than the actual number of nodes of that element. For one-dimensional problems, we 
simply interpolate the numerical solution onto 100 grid points within each finite element. These 
grid points are then utilized explicitly to visualize the finite element solution. For two-dimensional 



problems, equation (19) is again employed to evaluate the numerical solution at 256 unequally 
spaced grid points within each element. These points are the effective nodes associated with a 
spectral//ip finite element using a 15th order polynomial expansion (as geometrically characterized 
using the polynomial expansion associated with the actual element used in the numerical analysis). 
The refined set of points is then utilized to create a low-order visualization mesh that can be readily 
imported into Tecplot [I]. Alternatively, one can use any other procedure available in the literature 
for visualizing (scalar, vector and tensor) quantities interpolated using high-order approximations 
(for example, see [33j EZ] and references therein). 

We now illustrate the performance of p- and /i-refinements on various one- and two-dimensional 
test problems. For two-dimensional problems, the regions of the violation of the non-negative 
constraint are indicated using white color. 

5.2. One-dimensional problem with zero forcing function. This test problem is taken from 
Reference |22) . which addressed the low-order approximation (i.e., p = 1). For completeness, we 
shall outline the problem. The domain is taken as Q := (0, 1) with zero forcing function and non- 
zero Dirichlet boundary conditions. Mathematically, the test problem can be written as follows: 

d 2 c 

ac - -j-p = in (0, 1) (24a) 
c(x = 0) = c(x = 1) = 1 (24b) 

where a is a non-negative constant. The analytical solution in terms of a is given by 

1 - exp[-^a] r i i , exp[Va] - 1 \ r~ \ m\ 

c(x) = ri = . = exp^ax] + — expf-^ax] (25) 

exp[ 1 /aj — exp[— yaj exp|y ctj — exp|— yaj 



Herein we shall take a(x) = 1000. It should be noted that the solution becomes steeper near the 
boundary as a becomes larger but the solution is still infinitely differentiable. The computational 
mesh is obtained by discretizing the domain using four equal-sized finite elements. The concen- 
tration profiles using the single-field formulation and the two least-squares-based formulations are, 
respectively, shown in Figures [3j [4] and [5] The flux profiles obtained using the LS2 formulation 
(which is a least-squares-based mixed formulation) is shown in Figure Ril In Reference [22], it has 
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been shown that, for this test problem, the violation of the non-negative constraint vanishes with 
/i-refinement under the classical single-field formulation. From the numerical results presented 
in this subsection, one can conclude that the violation of the non-negative constraint also vanishes 
with p '-refinement under the classical single-field and least- squares formulations for one- dimensional 
problems with zero forcing function. 

5.3. One-dimensional problem with non-zero forcing function. This test problem is taken 
from Reference [31] . Consider pure diffusion (that is, decay is neglected) in domain := (—1,1) 
with homogeneous Dirichlet boundary conditions. The forcing function is taken to be 

/(x) = 200 exp[-10(x + 1)] (26) 

The test problem takes the following form: 

~S =/(x) inn (27a) 

c(x= -1) = c(x= 1) = (27b) 

The analytical solution is given by 

c(x) = -2 exp[-10(x + 1)] - (1 - exp[-20])x + (1 + exp[-20]) (28) 

We shall mesh the computational domain using one element, and solve the problem for various 
polynomial approximations in the element. For this test problem, the LSI and LS2 formulations 
are identical as a(x) = and D(x) = I. The obtained concentration profiles for various polynomial 
approximations under the classical single-field formulation and the least-squares formulations are 
shown in Figure [7j Although for this problem /(x) > Vx 6 Q, the classical single- field formulation 
(for p = 3) and the least-squares formulation (for p = 3 and p = 5) violate the non-negative 
constraint. The reason for the violation is that, under the classical single- field formulation, the Li 
projection of the forcing function onto the space spanned by p = 3 has negative values near x = 1. 
A similar reason holds for the least-squares formulation. 

5.4. Two-dimensional isotropic diffusion. This test problem was proposed by Burman and 
Ern [6]. The computational domain is a rectangle f2 := (0, 1) x (0,0.3) with homogeneous Dirichlet 
boundary conditions. The decay coefficient a(x) is taken to be zero. The diffusivity tensor is the 
second-order identity tensor. That is, 

D(x) = I (29) 

The forcing function is taken as follows: 

f 1 if (x, y) £ [0,0.5] x [0,0.075] , 
/(x) = \ J (30) 

I otherwise 
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Table 1. Two-dimensional isotropic diffusion: This table shows the variation of 
the minimum concentration with respect to p-refinement and /i-refinement. Since 
the medium is isotropic, the violation of the non- negative constraint vanishes with 
/i-refinement. Although the forcing function is non- negative, the violation of the 
non- negative constraint did not vanish with p-refinement. One of the reasons is that 
the L2 projection of the forcing function onto the polynomial space is not non- 
negative. 



h/p order 


p-refinement 


/i-refinement 


Single-field 


Least-squares 


Single-field 


Least-squares 


1 














2 


-1.35 x 10~ 4 


-1.41 x 10~ 4 


-8.84 x 10~ 5 


-1.08 x 10" 4 


3 


-3.13 x 1(T 5 


-3.17 x 10~ 5 


-1.81 x 10~ 5 


-2.57 x 10~ 5 


4 


-2.04 x 10~ 5 


-2.08 x 10~ 5 








5 


-4.42 x 10~ 6 


-4.46 x 10~ 6 








6 


-2.77 x 10~ 6 


-2.79 x 10~ 6 








7 


-1.99 x 10~ 7 


-2.01 x 10~ 7 








8 


-6.22 x 10~ 7 


-6.24 x 10~ 7 








9 


-2.60 x 10~ 8 


-2.60 x 10~ 8 








10 


-2.06 x 10~ 7 


-2.06 x 10~ 7 









Since the forcing function /(x) > in $7 and c p (x) = on d£l, from the maximum principle given 
in Theorem [l] we have c(x) > in Cl. It should be noted that since the diffusivity tensor is isotropic 



and there is no decay, both the least-squares formulations are identical (see equation (13)). The 
computational and visualization meshes are shown in Figure [8j The concentration profiles obtained 
under the single-field and the least-squares formulations for various polynomial approximations are, 
respectively, shown in Figures [9] and 10 



The minimum concentration under the single-field and the least-squares formulations for various 
order of p- and /i-refinements are given in Table [TJ As one can see, there is violation of the non- 
negative constraint under p-refinement, and the violation vanishes under /i-refinement. As discussed 
towards the end of Section [2j the reason is that the L2 projection of the forcing function on the 
polynomial space need not be non-negative under high-order approximations (even though the 
forcing function is non- negative) . Note that the L2 projection of the forcing function onto p = 1 
polynomial space is non-negative. 

5.5. Non-uniform anisotropic media. This test problem was originally proposed by Le Potier 
, and has been employed in many other numerical studies on maximum principles and the 
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non- negative constraint using the low-order approximation (e.g., see References [2H 123]). The test 
problem assumes a(x) = 0. The computational domain is a square Q := (0,0.5) x (0,0.5) with 
homogeneous Dirichlet boundary conditions enforced on the whole boundary. The diffusivity tensor 
is given by 

I y 2 + ex 2 -(1 - e)xy \ , 

d x = ; \ 7 (31) 

V -(l-e)xy x 2 + ey 2 / 

with e = 10 -3 . The forcing function is given by 

fl (x, y) G [0.125,0.375] x [0.125,0.375] , 
/(x) = <^ J L J (32) 

^ otherwise 

Two different meshes (one is a structured mesh and the other is an unstructured mesh) are employed, 



and are shown in Figures 11 and 12 The performance of the single-field, LSI and LS2 formulations 



using the structured mesh under both p- and /i-refinements is illustrated in Figures 13, 14 and 15 



The performance of these formulations using the the unstructured mesh is illustrated in Figures 



16 17 and 18 The minimum concentration under these three formulations is shown in Figure 
19 Since the anisotropy is strong, there will be violation of the non-negative constraint even on 
structured mesh for both p- and /i-refinements. For this test problem, it should be noted that LSI 
formulation did not perform as well as LS2 formulation. This is expected for problems involving 
heterogeneous and anisotropic media as the weight (s) in defining least-squares functional play a 



crucial role in the accuracy of the numerical results. (See equation ( 13 ) to note the different weights 
used in constructing LSI and LS2 formulations.) 

5.6. Anisotropic diffusion in a square domain with a hole. This problem has been used 
in References |21] [T9l [23j [22] with respect to the enforcement of the non-negative constraint but 
in the context of low-order approximation. The computational domain is a bi-unit square with a 
square hole of dimension [4/9,5/9] x [4/9,5/9]. The forcing function is taken to be /(x) = 0, and 
the decay coefficient is a(x) = 0. On the external boundary c p (x) = is prescribed, and on the 
internal boundary c p (x) = 2 is prescribed. The diffusivity tensor is given by 

D(x)= ( C ° S{9) S ' m(e) )( h ° )( COS{0) ~ SinW ] (33) 



sin(0) cos(0) j V k 2 J V sin(0) cos(0) 



Herein, we have taken = tt/6, and considered two different sets for diffusivity coefficients: 
(feijfo) = (1,100) and (fci,/c2) = (1,10000). The computational and visualization meshes are 
shown in Figure [20j The concentration profiles under the single-field, LSI and LS2 formulations 



for various polynomial approximations are, respectively, shown in Figures 21, 22 and 23 Figure 24 
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shows the minimum concentration under these three formulations for both sets of diffusivity coeffi- 
cients. From this figure, it is evident that the greater is disparity between the diffusivity coefficients 
k\ and k^ the greater is the violation of the non-negative constraint. Moreover, the extent of the 
violation did not decrease with with p- and /i-refinements for strong anisotropic medium (see the 
case ki = 1 and k<± = 10000). Hence, one can conclude that for strong anisotropic medium, p- and 
/i-refinement do not eliminate the violation of the non- negative constraint. 



6. CONCLUSIONS 

In this paper, we considered the diffusion-type equation, which is a second-order elliptic partial 
differential equation. This particular equation is known to satisfy maximum principles and the 
non-negative constraint under certain conditions. Many popular numerical formulations do not 
satisfy either maximum principles or the non- negative constraint. In the literature, it has been 
documented the performance of finite element formulations with respect to maximum principles for 
low-order elements. Herein, we considered the classical single- field formulation as well as two least- 
squares-based mixed formulations. We have systematically documented the performance of these 
formulations with respect to maximum principles and the non-negative constraint under high-order 
approximations. The main findings of this paper can be summarized as follows: 

(a) For one-dimensional problems, it is well-known that a uniform mesh is sufficient to satisfy the 
non-negative constraint and maximum principles under the low-order approximation. (For non- 
zero decay coefficient, the size of the element should be smaller than a critical size.) However, 



a uniform mesh is not sufficient for high-order approximations. (See subsections 5.2 and 5.3 



(b) For isotropic two-dimensional problems, a well-centered triangular mesh or a mesh with rectan- 
gular elements with aspect ratio between 1 / y/2 and \f2 is sufficient for the low-order approxi- 
mation to satisfy the non-negative constraint and maximum principles. Again, these conditions 



are not sufficient for high-order approximations. (See subsection 5.4 ) 
(c) For anisotropic two-dimensional problems, any restrictions on the mesh will not be sufficient, 
in general, to satisfy the non-negative constraint and maximum principles for both low-order 



and high-order approximations. (See subsections 5.5 and 5.6.) 
(d) The performance of a least-squares formulation depends on the weight (s) used in defining the 



least-squares functional. (See subsection 5.5, and Figures 14 and 18 



We shall conclude this paper with the following statement with a hope that it will motivate applied 
mathematicians, computational mechanicians, and numerical analysts to work on an interesting 
problem: A finite- element-based formulation or methodology that satisfies the non-negative con- 
straint and maximum principles for anisotropic diffusion on general computational grids under 
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high-order approximations is currently an unsolved problem with many important applications in 
engineering and applied sciences. 
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Figure 1 . High polynomial order one-dimensional C° Lagrange interpolation func- 
tions (cases shown are for p = 6) with (a) equal spacing of the element nodes, and 
(b) unequal nodal spacing associated with Gauss-Lobatto-Legendre points. 
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Figure 2 . Examples of various high-order two-dimensional master elements Vt e used 
in the present formulation: (a) four- node element (p = 1), (b) nine-node element 
(p = 2), (c) twenty-five- node element (p = 4), and (d) eighty-one-node element 
= 8). 
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-0.4 r 1 .0.4 r 1 

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 

(a) x (b) x 

Figure 3. One-dimensional problem with zero forcing function: The figure illus- 
trates the performance of the classical single-field formulation (which is also referred 
to as the Galerkin formulation) with respect to the non- negative constraint. In 
the right figure, we show the classical single-field formulation with respect to p- 
refinement. For comparison, the left figure shows the performance for an equivalent 
/i-refinement. For this test problem, the violation of the non-negative constraint 
decreases with both h- and p-refinements under the classical single-field formulation. 




-0.4 r i -0.4 r i 

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 



(a) * (b) 

Figure 4. One-dimensional problem with zero forcing function: The figure illus- 
trates the performance of the first-form of the least-squares formulation (LSI formu- 
lation) with respect to the non-negative constraint. In the right figure, we show the 
performance of the formulation with respect to p-refinement. For comparison, the 
left figure shows the performance for an equivalent /i-refinement. For this problem, 
the violation of the non- negative constraint decreases with both h- and p-refinements 
under the LSI formulation. 
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Figure 5. One-dimensional problem with zero forcing function: The figure illus- 
trates the performance of the second-form of the least-squares formulation (LS2 
formulation). In the right figure, we show the performance of the formulation due 
to p-refinement. For comparison, the left figure shows the performance for equiv- 
alent /i-refinement. For this problem, the violation of the non-negative constraint 
decreases with both h- and p-refinements under the LSI formulation. 




(a) x (b) x 

Figure 6. One-dimensional problem with zero forcing function: The figure shows 
the flux obtained using the second-form of the least-squares formulation (LS2 for- 
mulation). In the right figure, we show the performance of the formulation due to 
p-refinement. For comparison, the left figure shows the performance for equivalent 
/i-refinement. 
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Figure 7. One-dimensional problem with non-zero forcing function: This figure 
shows the obtained concentration profiles for various orders of polynomial approxi- 
mations under the classical single-field formulation (top figure) and the least-squares 
formulation (bottom figure). Although /(x) > Vx 6 both the classical single- 
field formulation (for p = 3) and the least-squares formulation (for p = 3 and p = 5) 
violate the non- negative constraint. 
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Figure 8. Two-dimensional isotropic diffusion: The top figure shows the mesh 
used in the numerical simulation, and the bottom figure shows the corresponding 
mesh used for visualization under high-order approximations. For more details on 



visualizing results using high-order approximations see subsection 5.1 
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Figure 9. Two-dimensional isotropic diffusion: This figure shows the contours of 
the concentration obtained using the single-field formulation for various polynomial 
approximations. The top figures are for p = 1 and p = 2, and the bottom figures are 
for p = 5 and p = 10. The regions in which the non-negative constraint is violated 
is shown in white color. 
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Figure 10. Two-dimensional isotropic diffusion: This figure shows the contours of 
the concentration obtained using the least-squares formulation for various polyno- 
mial approximations. The top figures are for p = 1 and p = 2, and the bottom 
figures are for p = 5 and p = 10. (Since the diffusivity is isotropic and decay coef- 
ficient is zero, both LSI and LS2 formulations are identical.) The regions in which 
the non-negative constraint is violated is shown in white color. 
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FIGURE 11. Non- uniform anisotropic media: The top figure shows the structured 
mesh used in the numerical simulation, and the bottom figure shows the correspond- 
ing mesh used for visualization under high-order approximations. For more details 



on visualizing results using high-order approximations see subsection 5.1 
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Figure 12. Non-uniform anisotropic media: The top figure shows the unstructured 
mesh used in the numerical simulation, and the bottom figure shows the correspond- 
ing mesh used for visualization under high-order approximations. For more details 



on visualizing results using high-order approximations see subsection 5.1 
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Figure 13. Non- uniform anisotropic media: This figure shows the concentration 
contours obtained using the single-field formulation and structured mesh (which is 



show in Figure 11). The left figures are using /i-refinement, and the right figures 



are using p-refinement. The top figures are for h/p = 2, the middle figures are for 
h/p = 5, and the bottom two figures are&r h/p = 10. 
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Figure 14. Non- uniform anisotropic media: This figure shows the concentration 
contours obtained using the LSI formulation and structured mesh (which is shown 



in Figure 11 ). The left figures are using h- refinement, and the right figures are using 



p-refinement. The top figures are for h/p = 2, the middle figures are for h/p = 5, 
and the bottom two figures are for h/p =§40. 
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Figure 15. Non- uniform anisotropic media: This figure shows the concentration 
contours obtained using the LS2 formulation and structured mesh (which is shown 



in Figure 11 ). The left figures are using h- refinement, and the right figures are using 



p-refinement. The top figures are for h/p = 2, the middle figures are for h/p = 5, 
and the bottom two figures are for h/p =§510. 
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Figure 16. Non- uniform anisotropic media: This figure shows the concentration 
contours obtained using the single-field formulation and unstructured mesh (which 



is shown in Figure 12). The left figures are using /i-refinement, and the right figures 



are using p-refinement. The top figures are for h/p = 2, the middle figures are for 
h/p = 5, and the bottom two figures areifcr h/p = 10. 
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Figure 17. Non- uniform anisotropic media: This figure shows the concentration 
contours obtained using the LSI formulation and unstructured mesh (which is shown 



in Figure 12). The left figures are using h- refinement, and the right figures are using 



p-refinement. The top figures are for h/p = 2, the middle figures are for h/p = 5, 
and the bottom two figures are for h/p =§7l0. 
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Figure 18. Non- uniform anisotropic media: This figure shows the concentration 
contours obtained using the LS2 formulation and unstructured mesh (which is shown 



in Figure 12). The left figures are using h- refinement, and the right figures are using 



p-refinement. The top figures are for h/p = 2, the middle figures are for h/p = 5, 
and the bottom two figures are for h/p =§810. 
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Figure 19. Non- uniform anisotropic media: This figure shows the minimum con- 
centration for various order of polynomial approximation for single-field, LSI and 
LS2 formulations using structured (top figure) and unstructured (bottom figure) 
meshes. According to the maximum principle, the concentration should be non- 
negative and should occur on the boundary. However, for all the three formulations 
the obtained minimum concentration is negative. The violation of the non-negative 
constraint did not vanish under either p- 9 or h- refinement. For better comparison, 
we have taken the y-axis to be the logarithm of the negative of the minimum con- 
centration. 
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Figure 20. Anisotropic diffusion in a square domain with a hole: The top figure 
shows the mesh used in the numerical simulation, and the bottom figure shows the 
corresponding mesh used for visualization under high-order approximations. For 
more details on visualizing results using high-order approximations see subsection 
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Figure 21. Anisotropic diffusion in a square domain with a hole: This figure shows 
the concentration contours obtained using the single-field formulation for k% = 1 
and = 10000. The left figures are for h- refinement, and the right figures are 
p-refinement. The top figures are for h/p = 2, the middle figures are for h/p = 5, 
and the bottom two figures are for h/p =4il0. 
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Figure 22. Anisotropic diffusion in a square domain with a hole: This figure shows 
the concentration contours obtained using the LSI formulation for k\ = 1 and 
&2 = 10000. The left figures are using /i-refinement, and the right figures using 
p-refinement. The top figures are for h/p = 2, the middle figures are for h/p = 5, 
and the bottom two figures are for h/p =42L0. 
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Figure 23. Anisotropic diffusion in a square domain with a hole: This figure shows 
the concentration contours obtained using the LS2 formulation for k\ = 1 and 
&2 = 10000. The left figures are for /i-refinement, and the right figures are p- 
refinement. The top figures are for h/p = 2, the middle figures are for h/p = 5, and 
the bottom two figures are for h/p = 10.43 
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Figure 24. Anisotropic diffusion in a square domain with a hole: This figure shows 
the minimum concentration under single-field, LSI and LS2 formulations for two dif- 
ferent cases: k\ = 1 and k<i = 100 (top figure), and k\ = 1 and ki = 10000 (bottom 
figure). We have shown the minimum concentration under various levels of p- and 
/i-refinements. For all the three formulations and under various levels of p- and 
/i-refinements, the minimum concentration is negative. For better comparison, we 
have taken the y-axis to be the logarithrn of the negative of the minimum concen- 
tration. (Note that the minimum concentrations obtain using the single-field and 
LS2 formulations are approximately the same under both p- and h- refinements.) 



